function h0 = calcH0_VMD(params)


h0 = zeros(params.mesh.ne,params.ntx);
for ii = 1:params.ntx
    
    rxXf = params.mesh.Xf(:,1) - params.loc(ii,1);
    ryXf = params.mesh.Xf(:,2) - params.loc(ii,2);
    rzXf = params.mesh.Xf(:,3) - params.loc(ii,3);
    rXf = sqrt(rxXf.^2 + ryXf.^2 + rzXf.^2);
    Ax = -ryXf./(4*pi*rXf.^3);
    
    rxYf = params.mesh.Yf(:,1) - params.loc(ii,1);
    ryYf = params.mesh.Yf(:,2) - params.loc(ii,2);
    rzYf = params.mesh.Yf(:,3) - params.loc(ii,3);
    rYf = sqrt(rxYf.^2 + ryYf.^2 + rzYf.^2);
    Ay = rxYf./(4*pi*rYf.^3);
    
    Az = zeros(size(params.mesh.Zf,1),1);
    
    A = [Ax;Ay;Az];
    h0(:,ii) = params.mesh.CURL'*A;
end


